{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "37a622fc-6c7a-4b52-909e-b6bcc5e0b2e5",
   "metadata": {},
   "source": [
    "# Quantum-Selected Configuration Interaction (QSCI)\n",
    "\n",
    "Quantum-Selected Configuration Interaction (QSCI) [1] is a quantum-classical hybrid algorithm designed for electronic structure calculations, particularly to compoute ground state and excited state energies of molecular systems. This method is one of the quantum subspace diagonalization (QSD). It is well-suited for NISQ-era quantum devices, focusing on using quantum resources efficiently while maintaining high accuracy.\n",
    "\n",
    "![quantum-subspace-diagonalization](images/qsd.png)\n",
    "\n",
    "QSCI selects this vector space from a sampled computational basis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1388e9b0-1b9b-4df4-890c-117622b771c4",
   "metadata": {},
   "source": [
    "## 0. Problem definition\n",
    "\n",
    "Let's define Hamiltonian of $\\mathrm{H_4}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "175153b6-3bb4-4378-a174-cead14589e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --quiet --break-system-packages openfermionpyscf==0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "93c9a86f-6eab-48fb-9edf-51c7d25e9f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openfermion\n",
    "import openfermionpyscf\n",
    "from openfermion.transforms import jordan_wigner, get_fermion_operator\n",
    "import cudaq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "336fd99f-7e4e-4401-a60e-eecb9c869299",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of hydrogen atoms.\n",
    "hydrogen_count = 4\n",
    "\n",
    "# Distance between the atoms in Angstroms.\n",
    "bond_distance = 0.7474\n",
    "\n",
    "# Define a linear chain of Hydrogen atoms\n",
    "geometry = [(\"H\", (0, 0, i * bond_distance)) for i in range(hydrogen_count)]\n",
    "\n",
    "basis = \"sto3g\"\n",
    "multiplicity = 1\n",
    "charge = 0\n",
    "\n",
    "molecule = openfermionpyscf.run_pyscf(\n",
    "    openfermion.MolecularData(geometry, basis, multiplicity, charge), run_fci=True\n",
    ")\n",
    "molecular_hamiltonian = molecule.get_molecular_hamiltonian()\n",
    "fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian)\n",
    "qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)\n",
    "qubit_hamiltonian.compress()\n",
    "\n",
    "hamiltonian = cudaq.SpinOperator(qubit_hamiltonian)\n",
    "\n",
    "num_qubits = hamiltonian.qubit_count\n",
    "electron_count = molecule.n_electrons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "908956ec-a201-416d-bc37-f09c6ba02b85",
   "metadata": {},
   "source": [
    "## 1. Prepare an Approximate Quantum State\n",
    "\n",
    "Start with a variational ansatz (e.g., from VQE) to approximate the ground state wavefunction on a quantum computer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "04b6e475-72bc-4f4a-bed7-32488d044e0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def kernel(thetas: list[float]):\n",
    "\n",
    "    qubits = cudaq.qvector(num_qubits)\n",
    "\n",
    "    # Hartree-Fock\n",
    "    for i in range(electron_count):\n",
    "        x(qubits[i])\n",
    "    # UCCSD\n",
    "    cudaq.kernels.uccsd(qubits, thetas, electron_count, num_qubits)\n",
    "\n",
    "\n",
    "parameter_count = cudaq.kernels.uccsd_num_parameters(electron_count, num_qubits)\n",
    "\n",
    "\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "\n",
    "# Define a function to minimize\n",
    "def cost(thetas: list[float]):\n",
    "    exp_val = cudaq.observe(kernel, hamiltonian, thetas).expectation()\n",
    "    return exp_val\n",
    "\n",
    "\n",
    "exp_vals = []\n",
    "\n",
    "\n",
    "def callback(xk):\n",
    "    exp_vals.append(cost(xk))\n",
    "\n",
    "\n",
    "# Initial variational parameters.\n",
    "np.random.seed(15)\n",
    "x0 = np.random.normal(0, np.pi, parameter_count)\n",
    "\n",
    "# Use the scipy optimizer to minimize the function of interest\n",
    "vqe_result = minimize(cost, x0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02dc831a-c41f-4abf-8347-4df7a0d5ca72",
   "metadata": {},
   "source": [
    "## 2 Quantum Sampling to Select Configuration\n",
    "\n",
    "Measure this quantum state multiple times to obtain bitstrings, which correspond to Slater determinants (electronic configurations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0eb10a37-9384-41f5-bb9e-92d91700434a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ 00001111:2 00011011:12 00011110:18 00100111:27 00101101:25 00110011:2 00110110:22 00111001:193 01011010:13 01101100:23 01110010:2 10000111:34 10001101:115 10010011:18 10010110:5 10011100:11 10100101:40 10110001:4 11000011:170 11001001:16 11001100:8 11010010:46 11011000:137 11100001:7 11100100:28 11110000:22 }\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sample_result = cudaq.sample(kernel, vqe_result.x)\n",
    "print(sample_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "836fb668-4bb7-4214-aaaf-72d8f697a384",
   "metadata": {},
   "source": [
    "## 3. Classical Diagonalization on the Selected Subspace\n",
    "\n",
    "Construct a truncated Hamiltonian matrix using only the selected determinants.\n",
    "Diagonalize this matrix classically to obtain improved energy estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "eb220374-44ca-4183-9e1a-f27cc62e8fd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_element(\n",
    "    operator: cudaq.SpinOperator, b1: str | np.ndarray, b2: str | np.ndarray\n",
    ") -> complex:\n",
    "    \"\"\"\n",
    "    Compute an inner product <b1|operator|b2>.\n",
    "    \"\"\"\n",
    "\n",
    "    if isinstance(b1, str):\n",
    "        b1 = np.array([i == \"1\" for i in b1], dtype=bool)\n",
    "    if isinstance(b2, str):\n",
    "        b2 = np.array([i == \"1\" for i in b2], dtype=bool)\n",
    "\n",
    "    num_qubits = hamiltonian.qubit_count\n",
    "    num_terms = hamiltonian.term_count\n",
    "\n",
    "    bsv = np.zeros((num_terms, 2 * num_qubits), dtype=bool)\n",
    "    coeffs = np.zeros(num_terms, dtype=np.complex128)\n",
    "    for i, term in enumerate(operator):\n",
    "        coeffs[i] = term.evaluate_coefficient()\n",
    "        bsf = term.get_binary_symplectic_form()\n",
    "        n = len(bsf) // 2\n",
    "        for j, e in enumerate(bsf):\n",
    "            if j < n:\n",
    "                bsv[i, j] = e\n",
    "            else:\n",
    "                bsv[i, j - n + num_qubits] = e\n",
    "    x = bsv[:, :num_qubits]\n",
    "    z = bsv[:, num_qubits:]\n",
    "\n",
    "    # b1 ⊕ b2 must match x for nonzero contribution\n",
    "    delta = np.bitwise_xor(b1, b2)\n",
    "    match = np.all(x == delta, axis=1)\n",
    "\n",
    "    if not np.any(match):\n",
    "        return 0.0\n",
    "\n",
    "    matched_x = x[match]\n",
    "    matched_z = z[match]\n",
    "    matched_coeffs = coeffs[match]\n",
    "\n",
    "    # i^{x·z}\n",
    "    x_dot_z = np.sum(np.bitwise_and(matched_x, matched_z), axis=1) % 4\n",
    "    i_phases = np.array([1, 1j, -1, -1j], dtype=np.complex128)[x_dot_z]\n",
    "\n",
    "    # (-1)^{b1·z}\n",
    "    b1_dot_z_parity = np.sum(np.bitwise_and(b1, matched_z), axis=1) % 2\n",
    "    sign_phases = np.where(b1_dot_z_parity == 0, 1, -1)\n",
    "\n",
    "    result_terms = matched_coeffs * i_phases * sign_phases\n",
    "    return np.sum(result_terms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5c095c67-1222-4366-bf0b-4a4d847a2841",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "\n",
    "bitstrs = list(k for k, _ in sample_result.items())\n",
    "subspace_dimension = len(bitstrs)\n",
    "\n",
    "values = []\n",
    "row_ids = []\n",
    "column_ids = []\n",
    "\n",
    "for i in range(subspace_dimension):\n",
    "    for j in range(i, subspace_dimension):\n",
    "        elem = pauli_element(hamiltonian, bitstrs[i], bitstrs[j])\n",
    "        if abs(elem) > 1e-6:\n",
    "            values.append(elem)\n",
    "            row_ids.append(i)\n",
    "            column_ids.append(j)\n",
    "            if i != j:\n",
    "                values.append(elem.conjugate())\n",
    "                row_ids.append(j)\n",
    "                column_ids.append(i)\n",
    "\n",
    "values = np.real_if_close(values)\n",
    "\n",
    "truncated_hamiltonian = scipy.sparse.coo_array(\n",
    "    (values, (row_ids, column_ids)), shape=(subspace_dimension, subspace_dimension)\n",
    ")\n",
    "eigvals, _ = scipy.sparse.linalg.eigsh(truncated_hamiltonian, 1, which=\"SA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb40b6c3-063a-45df-8867-747fac2bb354",
   "metadata": {},
   "source": [
    "## 5. Compuare results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "827cfd92-9208-4057-959f-db70853a32d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.1318697852134716"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## QSCI\n",
    "eigvals[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "df643aee-1967-4d3d-aebd-dda13892c9ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.143554580446029"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## FCI Energy\n",
    "molecule.fci_energy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd08b5ed-dc5a-4f25-86fd-a6c805ee8ee6",
   "metadata": {},
   "source": [
    "## Reference\n",
    "\n",
    "[1] \"Quantum-selected configuration interaction: classical diagonalization of hamiltonians in subspaces selected by quantum computers\" K. Kanno, M. Kohda, R. Imai, S. Koh, K. Mitarai, W. Mizukami, and Y. O. Nakagawa (2023) arXiv: 2302.11320"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
